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Abstract 

The simplest solvable problem of stress transmission through a static granular ma- 
terial is when the grains are perfectly rigid and have an average coordination number 
of z — d+l. Under these conditions there exists an analysis of stress which is indepen- 
dent of the analysis of strain and the d equations of force balance Vj dij (r) = gi (r) 
have to be supported by equations. These equations are of purely geometric 

origin. A method of deriving them has been proposed in an earlier paper [|| . In this 
paper alternative derivations are discussed and the problem of the " missing equations" 
is posed as a geometrical puzzle which has yet to find a systematic solution as against 
sensible but fundamentally arbitrary approaches. 

1 Introduction 

Granular media are ubiquitous yet complex materials with puzzling properties The 
simplest model of a static granular material is that where grains are considered to be 
perfectly hard, perefectly rough and each grain a has a coordination number = d + 1 
(where d is the dimension of the system). Under these conditions Newton's equations 
of intergranular force and couple balance can be solved IQ . The system is in the state of 
mechanical equihbrium and particles can not experience deformation under load so there is 
no displacement field present. Thus the only immediately available macroscopic equation 
has the form 



ViCTij(r) = gi{r), (1) 

where o"jj(r) is the macroscopic stress tensor and gi{r) is external force at the boundaries. 
The vector equation ([l|) gives d equations for '^^'^^^^ components of o"jj(r) leaving 
further equations required to solve for the macroscopic stress tensor. Thus one should 
be able to derive them from the geometry of the contact network which is assumed to 
be specified. In order to put the above comments into formulae we draw a diagram (see 
Figure |l|) of one grain in contact with 4 nearest neighbours (i.e. grains that are in contact 
with the reference grain a). 

The geometrical specification of the system is given by the set of — contact points 
{R"^*}. The centroid of contacts of the reference grain a is defined by vector R" 



R' 



Figure 1: Cross-section of the first coordination shell of the reference particle a 



where the summation sign means the sum over all nearest neighbours of the grain a. 
The distance between grains a and P is defined as the distance between their centroids of 
contacts 

R"/5 = R/3 _ R° = - r^", (3) 

where r"^ is the vector joining the centroid of contact with the contact point. The second 
vector which characterises the relative position of neighbouring centroid with respect to 
the contact point is defined by (see Figure ^) 




Figure 2: Segment of the contact network between nearest neighbours a and /? 



Newton's laws of force and couple balance for every grain give us the system of 
equations for interparticle forces f "'^ (see Figure ^) 







(5) 



/r + /r" = , (6) 

Y.'^Ulffr-^=cr. (7) 
where is the external force acting on grain a and c° is the external couple. The 




Figure 3: Intergranular forces and the local geometry of two grains in contact 
tensorial force moment for grain a is defined as 

and it is symmetric tensor given that c" = 0. The macroscopic stress tensor Uij{r) can 
be obtained by averaging over the packing 

a^.ir) = { St,) , (9) 

where {■■■) = ^ S^=i • ■ ■ '^(r — R") in the simplest case and V is the volume of the pack- 
ing. The method offered by the authors in was to consider the probability functional 
for the set {S^A 




where 



a,l3 (3 

^KY^^^Ulffr^^) (11) 



X 



and the normaUsations, N and M are functions of contact network configuration. The 
main goal of Ref . 0] was to transform (p!0|) into 



p{ss] = nn^^^S^5 - 9t)5[Y.prAs^i) (12) 

where delta-functions contain the complete system of equations for the set of tensorial 
force moments {5*^"}. The method of Ref.0] was to exponentiate all delta- functions in 

(0) 

N 

p{5^«}= / JJ e*^P/"^PC"2^7°^A"p7?°^, (13) 



where 



+rir[fr + ft 

; /°^, Af and r/f'^ 
for the conjugate fields {Cj?} and {7^"} was obtained 



a f3 

^ (14) 



After integrating out the fields f°'^ , Af and r/"'^ the following system of linear equations 



C."rf - 7f = C,?rf - 7f . (15) 



From (15) it was shown in Ref.[g] that has the representation 



Q = Q^ + Q*, (16) 
where a particular solution Q'j^ gave the first delta- function in ( [T^ ) 

and a complimentary function * which satisfies the following system of linear equations 



C.r^^'-Cfrf = 0. (18) 

This system of linear equations for {Cij*} gave the required -^'^(^~^) constraints on {S^}. 
The present paper concentrates on C° * and uses the result of Ref . [|| that 



TV 

P {S^'J} = n K ^ - - 5.") P {S^beomeiry} , (19) 

a=l /3 

where P ^S^\geometry^ contains the set of ^^'^2 missing equations which in continuum 
limit might take a form 



Pijkl^^kl + ^ jTijkl<ykl + "^j^lUijklCTkm + ••• — . (20) 

In the first term we have Pjjfc; = —Pjj/j/ which gives the correct number of equations. The 
missing equation ( pO|) is of purely geometric origin. We deliberately avoid using the term 
"constitutive relation" (for there is no deformation or displacement in this model) and call 
Eq. ( [20| ) the stress-geometry equation. The authors believe this kind of situation offeres 
a new kind of challenge in theoretical physics and although this paper presents a solution, 
it does not have the elegance and completeness that one might expect from the solution 
derived by means of some variational principle. There may be some analogy here to the 
problems of dynamics where, although the Lagrange equations with Lagrange multipliers 
will solve any non-holonomic problem, the really powerful way is to use the Gibbs-Appell 
equations |^]. 



2 The missing stress-geometry equation 

In order to obtain P ^S^^lgeometry^ and derive ( pO|) we need to solve Eq.(18). Let us 
describe the method of Ref.Q. We notice that in order to obtain the precise number of 
missing equations (which is '^'^^2 ^^ )) Eq. (|T^) appears to be too many equations. Because 
{r"^} satisfy the linear relation r"'^ = from Eq.(|^) there are many internal 
identities and careful counting shows that it can only contain Nd equations. For example 
if when solved Eq. (18) were to give 



/- a * 



5I2 



(21) 



then 



P {S,^\geometry} = / J] e^E^ (^fiCfl* +25i'^,Cir + 5,'|C.r) , (22) 

gives = and no constraint on . One can force Eq.(p!8|) into only two equations 
(for d = 2) by taking scalar product with vectors a°^^ and 6"^ and then summing over /3: 

E * - E * r/" = , (23) 

P /3 



f^af^ ^a* a/3 _ ,af3./3* /3. 

/3 13 



(24) 



which gives us 2N equations. Suppose that we regard the differences between Q^* and 



C,* to be expandable in R"^. 



yal3 



/-a* I p( 



+ 



(25) 



The first approximation is then 



and similarly from Eq.(|2 



/3 



al3 > 



ij 



a(3 



0, 



(26) 



/3 



(27) 



We have declared in Ref.[|| that there are two obvious vectors to be used in Eqs.(p6D and 
(|27|), namely R"^ and Q"^. Hence we obtain two configuration tensors (analogous, but 
different from the fabric tensors used in the soil mechanics literature 01) 



and 



/3 



af3 pa/3 



r -i OL 



(28) 



(29) 



then we have 



(^a* pa ^ Q 



^ij ij 



(30) 
(31) 



and after exponentiating Eqs.(30,31) in P {geometry^ and eliminating the auxilirary 
fields one finds the missing stress-geometry equation 



Q a 


^11 


a 

'^11 






Q a 
'-'12 


^12 


ri a 
^\2 


= 0. 


(32) 


Q a 
'-'22 


rp a 


/-I a 
'^22 







The first approximation ( p^j27| ) (christened the "first coordination shell approximation") 
can be illustrated in the following way: we rearrange Eq.([l^) and after multiplying by 
R"^ and Q"^ obtain 

Cr E ^r"" + - ^if ) = ' (33) 

13 13 

c.r E + E(c.r - cf ) rf Qf = o , (34) 

/3 /3 



It is clear that if the second term in Eqs.mS) and (|3J) is neglected one obtains Eqs.(|30D 
and (^Tj). A naive attempt to transform (|32D into the macroscopic equation for o"jj(r) by 
averaging S^, F-" and fails (in the case of an isotropic configuration) because 

(G.") = ^ E E Qt'Rf = E E = (35) 

a=l /3 13=1 a 



Thus for an isotropic packing the first term in Eq. ( |20| ) vanishes in the first coordination 
shell approximation. This gives rise to the conditional probability distribution functions. 
Thus if we are given S{2, the probability of finding Sf{ — S22 is 



and vice versa. This distribution can then be introduced for corresponding components 
of the macroscopic stress tensor subject to the absence of mesoscopic correlations in the 
packing. For an anisotropic configuration ( |32|) characterized by the distribution of {G^} 
with some nonvanishing average yields macroscopic equation when 



o"!! - (T22 = 2(Ti2tan(/). (37) 

where (f) is the angle of repose in the case when configuration is prepared in the form of 
a sandpile . The obvious criticism of the derivation method of Ref. [Q] is that one could 
employ some different vector to obtain Eqs.(^) and (|3l|). For example instead of R"^ 
one could use R"^k(R"'^, Q"'^) where k is any scalar function of R°^ and Q"^; or one 
could go to the next coordination shell of the reference grain a and employ R'^''' and Q^"^ 
(where 7's are the other two neighbours of /3, see Figure Q) and so on. Thus we have 
offered a path to the missing equation which works also in 3-D, but it is not unique 
Presumably the internal symmetries of (18) will lead to the same macroscopic equation 



when Eqs. (l26| , |27D are used for any a"'^ , 5"^, but it is not easy to see how. In the next 
section we will offer some new viewpoints which approach the problem from a different 
standpoint, but which will confirm the earlier results, and suggest new approaches. 



3 New methods of derivation 

In the previous section we have shown that it is possible to go from Newton's equations 
(|^]^ to equations for conjugate fields ( p!5| ) that can be then used to derive the complete 
set of equations for {S^^}. We can reverse the process and from 



r ^ 

P {St,\geometTy] = y [] e^^^ s^c^' 5{C,^*rf - C,f r/") VC* , (38) 



obtain 



P {St,\geometry] = / J] '^(^^? " E Pf'-f) ^^Pf + P^) ^P"^ • (39) 

Since the only constraint on {p"^} is given by 

pf+pP- = Q (40) 
the missing stress-geometry equation comes from the condition on S[j that 

st, = \i: {pf^f + pf^f) (41) 

can be satisfied by a set of " pseudo- forces" {p°^} which obey ([40|), but are free from the 
constraint of Newton's second law (^). " Pseudo- force" p"^^ can be written in the following 
form 



pf = + Rf + (x- - x^) Qf , (42) 

,nd provides Eq. (|3^). We can think of {S^} as a component 

vector S, and {p"'^} as a ^"^(^^^^ component vector P. Then the set of equations (pl| ) 
can be written as 



S = RP (43) 

where R is a x ^'^(^+^) matrix with a large repetition of elements since r = 0. 

This matrix has ^'^^'^ -"-^ zero eigenvalues which gives the right number of constraints on 
The argument here is of the "must be so" type. Apart from the easy proof that 
Det R = (by adding rows) we have not succeeded in proving that the ^'^(^+^) x ^'^(^+^) 
matrix has ^'^(^~^) zero eigenvalues, but it "must be so"! Another and perhaps easier 
method is to return to the Eq.(p^) for ^^"* and sum it over j3 using r"^ = 0. 
Therefore we have 



ECf-r = 0. (44) 



This can be used in Eq. (|39| 



N 



P {S,^j\geometry] = f f[ e'^^ ^S^^T CiJ ^C* , (45) 

as in the derivation of ([4l|) , but now losing a positional index 

/3 



As before this equation implies a relationship between the components of S^. Suppose 

now that we are looking for an ansatz for Since (jif is a vector it can be represented 
as a superposition of two obvious candidates R"^^ and Q"'^ 

where new quantities ip" and are scalars. This gives us 



where H^'j = Q°^^Q°^^. Eliminating -0" and x"^ leads to 



(47) 



(48) 
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(49) 
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Note that i/j" will on average be a multiple of F^" and this gives us Eq. (32). It is possible 
to obtain from ([46|)the third term in (pO|). The second term is non- vanishing only in special 
cases of periodic arrays P]. Let us construct the following interpolation of 4>l^ 



after substituiting it into ( ^6|) and summing it over /3 using r 

Vfc " 



a/3 



we have 



Q a 



(50) 



(51) 



Crude averaging gives us 

1 ^ 



(52) 



a p 



We can now eliminate (pi and obtain the well-known Navier equation which imposes kine- 
matic compatibility on the stresses 



+ 



yy 



dxdy 



0. 



(53) 



This equation corresponds to the third term in (pO|) in the case of an isotropic packing 
and implies that the stress tensor components can be expressed in terms of the Airy 
function whose discrete analogues are ip"- and in (47). This also means that the set 
of microscopic constraints (^) is consistent with the macroscopic equation (0) . However in 
3-D we have encountered mathematical difficulties with this simple averaging procedure. 
This highlights the puzzle we face. One could put much more complicated versions of 
<?!)j^, for example any functional form which employs vectors constructed out of the vectors 
of the Figure ^ or indeed the next coordination shell (see Figure ^). The problem is a 
geometric puzzle. A set of contact points {R°^*} corresponds to the packing of grains 
with coordination number z = d+\. From these points associated with the reference grain 




Figure 4: Propagation of coordination shell vectors through the contact network 



a one can construct the vectors of Figure |2|, and find the basic equations using Cij * or 
p"/^, or (j)^. With these equations one may take ad hoc steps to obtain the complete set of 
equations for the macroscopic stress tensor aij (r) , but we have failed to find a systematic 
procedure such as is possible for the corrections of the stress-force equation as outlined 
in Ref.Q. The problem is now purely geometric, but we emphasize that although real 
granular materials have many features omitted here, we are studying the simplest possible 
case and the geometric puzzle offered here although difficult is quite basic, and no simpler 
case can be found. 
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